print("#######################################################################")
print("# Title: Gene-level SAM                                               #")
print("# Author: Yongkee Cho (yongkeecho@wustl.edu)                          #")
print("# Date: Aug. 27 2012                                                  #")
print("# Description: This script runs probeset-level SAM analysis and       #")
print("#              returns gene-level SAM-t score (best among matching    #")
print("#              probesets                                              #")
print("# Requirements: impute, samr, bioC/hgu133plus2.db                     #")
print("#                                                                     #")
print("# Usage: Rscript gene_sam.r test.csv 'Two class unpaired' test_g.csv  #")
print("#                           test_fc.csv test_samt.csv                 #")
print("#######################################################################")

# reads arguments
args <- commandArgs(trailingOnly = TRUE)

# reads data set
print("Reading sam data file...")   
sam_csv <- as.matrix(read.csv(args[1], header = FALSE))
y <- as.integer(sam_csv[2, 2:ncol(sam_csv)])
ny <- length(y)
x <- matrix(as.double(sam_csv[3:nrow(sam_csv), 2:ncol(sam_csv)]), ncol = length(y), byrow = FALSE)
probesetid <- sam_csv[3:nrow(sam_csv), 1]

# executes SAM
print("Executing SAM...") 
if ("impute" %in% rownames(installed.packages()) == FALSE) {
   install.packages("impute")
}
if ("samr" %in% rownames(installed.packages()) == FALSE) {
   install.packages("samr")
}
samfit <- SAM(x, y, resp.type = args[2], geneid = probesetid, nperms = 100)
tt <- samfit$samr.obj$tt
foldchange <- samfit$samr.obj$foldchange

# maps probesetid into geneid
print("Mapping probeset id into gene id...")
if ("hgu133plus2.db" %in% rownames(installed.packages()) == FALSE) {
   source("http://bioconductor.org/biocLite.R")
   biocLite("hgu133plus2.db")
}
library(hgu133plus2.db)
map <- hgu133plus2ENTREZID
mappedprobes <- mappedkeys(map)
mappedlist <- as.list(map[mappedprobes])
nn <- length(mappedprobes)
idx <- rep(0, nn)
geneid <- rep("", nn)
for (i in 1 : nn) {
     idx[i] <- which(probesetid == mappedprobes[i])
     geneid[i] <- mappedlist[[i]]
}
ttmap <- tt[idx]
data <- data.frame(x = x[idx,], geneid = geneid, tt = ttmap)

# aggregate probeset-level ouput into gene-level output using minmax function
minmax <- function(x) {
  minx <- min(x)
  maxx <- max(x)
  if (abs(minx) > abs(maxx)) {
    return(minx);
  }
  else {
    return(maxx);
  }
}

aggr <- aggregate(tt ~ geneid, data = data, minmax)
ng <- length(aggr[,1])
idx <- rep(0, ng)
for (i in 1 : ng) {
  idx[i] <- which(data$tt == aggr[i,2])
}
geneidg <- aggr$geneid

# write gene-level expression value
xg <- x[idx,]
colnames(xg) <- y
rownames(xg) <- geneidg
write.table(xg, args[3], sep=",")

# write gene-level fold change value
foldchangeg <- log2(foldchange[idx])
write(paste(geneidg, foldchangeg, sep=","), args[4])

# write gene-level SAM t-statistics
ttg <- tt[idx]
write(paste(geneidg, ttg, sep=","), args[5])
